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Abstract 

Recent results from the Pierre Auger Observatory (PAO) indicate that the composition of ultra-high-energy cosmic rays (UHECRs) 
with energies above 10 19 eV may be dominated by heavy nuclei. An important question is whether the distribution of arrival 
directions for such UHECR nuclei can exhibit observable anisotropy or positional correlations with their astrophysical source 
objects despite the expected strong deflections by intervening magnetic fields. For this purpose, we have simulated the propagation 
of UHECR nuclei including models for both the extragalactic magnetic field (EGMF) and the Galactic magnetic field (GMF). 
We find that the GMF is particularly crucial for suppressing the anisotropy as well as source correlations. Assuming that only 
iron nuclei are injected steadily from sources with equal luminosity and spatially distributed according to the observed large scale 
structure in the local Universe, at the number of events published by the PAO so far (69 events above 5.5 x 10 19 eV), the arrival 
distribution of UHECRs would be consistent with no auto-correlation at 95% confidence if the mean number density of UHECR 
sources n s > 10~ 6 Mpc~ 3 , and consistent with no cross-correlation with sources within 95% errors for n s > 1(T 5 Mpc 3 . On 
the other hand, with 1000 events above 5.5 x 10 19 eV in the whole sky, next generation experiments can reveal auto-correlation 
with more than 99% probability even for n s < 10~ 3 Mpc -3 , and cross-correlation with sources with more than 99% probability for 
n s < 10 Mpc 3 . In addition, we find that the contribution of Centaurus A is required to reproduce the currently observed UHECR 
excess in the Centaurus region. Secondary protons generated by photodisintegration of primary heavy nuclei during propagation 
play a crucial role in all cases, and the resulting anisotropy at small angular scales should provide a strong hint of the source location 
if the maximum energies of the heavy nuclei are sufficiently high. 
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1. Introduction 

The origin of ultra-high-energy cosmic rays (UHECRs) with 
energies > 10 19 eV is an intriguing mystery in modern as- 
trophysics. Their sources are generally believed to be extra- 
galactic objects, although some Galactic objects may also be 
viable. In either case, the maximum energy of 10 20 eV can 
only be achieved in extreme environments [1]. Prominent 
source candidates suggested so far include active galactic nu- 
clei (AGN) [2, 3, 4, 5, 6, 7, 8, 9, 10], gamma-ray bursts (GRBs) 
[11, 12, 13, 14], neutron stars or magnetars [15, 16, 17, 18] 
and clusters of galaxies [19, 20, 21]. If UHECRs with en- 
ergies above ~ 6 x 10 19 eV are mainly protons, their propa- 
gation distance should be limited by interactions with cosmic 
microwave background (CMB) photons [22, 23, 24, 25, 26], 
the so-called Greisen-Zatsepin-Kuz'min (GZK) mechanism, so 
that they are observable only from sources sufficiently nearby 
(typically < 200 Mpc). Since all known astrophysical ob- 
jects are distributed inhomogeneously at such distances in the 
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local Universe, anisotropies in the distribution of UHECR ar- 
rival directions are expected as long as cosmic magnetic fields 
are weak enough to allow quasi-rectilinear propagation of the 
UHECRs. Searches for correlations between the arrival di- 
rections and the celestial positions of potential source candi- 
dates should then provide valuable clues to reveal their origin 
[27,28,29]. 

The Pierre Auger Observatory (PAO) operating in the south- 
ern hemisphere has reported possible evidence of a gradual in- 
crease in the average mass of UHECR particles with energies 
above 10 18 5 eV by analyzing the average depth of the shower 
maximum < X nYdx > and the root mean square of the shower-to- 
shower fluctuations of X max (RMS (X mdx )), assuming that cur- 
rent hadronic interaction models are realistic at these energies 
[30]. This is in contrast to measurements of similar quantities 
for the northern sky by the High Resolution Fly's Eye (HiRes) 
experiment that are consistent with an UHECR composition 
dominated by protons [31, 32]. If UHECRs consisted mainly 
of heavy nuclei rather than protons as indicated by the PAO re- 
sults, their deflections during propagation in magnetized extra- 
galactic and Galactic environments would be much larger than 
those for protons, and the hopes for observing anisotropy may 
be considerably weakened. On the other hand, PAO has seen 
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statistically significant correlations between the arrival direc- 
tions of UHECRs with energies above ~ 6 x 10 19 eV and the 
projected positions of AGN with z < 0.018 within an angular 
scale of 3.1° [33, 34] l . This does not necessarily point to AGN 
as the true sources, since they may only be tracers of the more 
numerous population of galaxies that make up the anisotropic, 
local cosmic matter distribution, and there are also consider- 
able uncertainties in the angular displacement due to deflec- 
tions by both the extragalactic magnetic field (EGMF) and the 
Galactic magnetic field (GMF). Nevertheless, this does imply 
that UHECR sources are associated with the matter distribu- 
tion. Independent analyses of the PAO data have confirmed the 
correlations of UHECR arrival directions with various kinds of 
nearby astrophysical objects [36, 37, 38, 39]. Such correlations 
were predicted in several theoretical studies under the assump- 
tion of protons as the primary UHECR particles [27, 28, 29]. 
We must beware that at the moment, the positional correla- 
tions (or anisotropy) are observed only above ~ 6 x 10 19 eV, 
whereas the PAO results on the composition are available only 
up to ~ 6 x 10 19 eV due to the lack of event statistics. We 
cannot rule out the possibility that protons become dominant 
above this energy and produce the anisotropy. However, there 
are presently no hints of such a sudden change in the compo- 
sition, and the simplest extrapolation of the measured trends of 
< ^max > and RMS(X miix ) to higher energies would entail a 
composition dominated by heavy nuclei. Since the appearance 
of anisotropies or cross correlations with astrophysical objects 
in the case of UHECR nuclei is not trivial, deeper, relevant stud- 
ies of UHECR propagation are warranted in order to achieve a 
consistent interpretation of the current data set as well as to 
quantify the prospects for future observations. 

Propagation of UHE nuclei has been studied previously by 
many authors from various perspectives [24, 25, 26, 40, 41, 42, 
43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 
60, 61]. Many have studied the photodisintegration process and 
consequent spectra in detail without accounting for magnetic 
fields, but only a few have focused on the effects of the magnetic 
fields. Ref. [44] calculated the propagation of UHE nuclei in 
a uniform turbulent EGMF model and examined their trajecto- 
ries and photodisintegration interactions in intergalactic photon 
fields, concluding that a characteristic feature in the spectrum 
may result. Refs. [47, 48] discussed the propagation of nuclei 
in a structured EGMF model obtained through a numerical sim- 
ulation of cosmological structure formation, showing that while 
heavy nuclei can be strongly deflected (> 20°) by the EGMF, 
nuclei (and their secondaries) with relatively small deflections 
can generate anisotropy on intermediate angular scales. 

In this paper, we study the propagation of UHECR nuclei in 
cosmic radiation and magnetic fields under different assump- 
tions for the source properties and the fields, focusing on the 
resulting distribution of arrival directions. In addition to a de- 



1 Contrastingly, HiRes has not found significant cross correlations with 
nearby extragalactic objects based on the same analysis method as PAO [35]. 
However, the two results cannot be naively compared because the exposure of 
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of the sky. 



tailed treatment of the photopair, photomeson and photodisin- 
tegration interactions of UHE nuclei with background photons, 
deflections in GMF and EGMF are taken into account for both 
primary nuclei as well as secondary nuclei arising from photo- 
disintegration. By statistically analyzing the calculated arrival 
distributions, we discuss the implications for current and future 
observations of anisotropy and cross correlations with sources. 

Motivated by the recent PAO result, we adopt a model for the 
UHECR source distribution that follows the density of large- 
scale structure actually observed in the local Universe. This is 
in contrast to Refs. [47, 48] who took the sources to be re- 
lated to the matter distribution in a numerical simulation, but 
not reflecting the actual Universe. Because photodisintegration 
limits the propagation distance of heavy nuclei in a way similar 
to the GZK mechanism for protons, a realistic account of local 
source inhomogeneity is essential. Unlike most previous work 
on the subject, here we also account for the GMF as well as the 
EGMF, since the former unavoidably affects all UHECRs ar- 
riving at Earth. Concerning the UHECR source composition, a 
wide variety of possibilities are currently allowed that are con- 
sistent with both the observed spectrum and the composition in- 
dicators X mux and RMS(X m . dx ) (e.g., [60]). This study assumes 
a pure iron composition, which can also reproduce the highest- 
energy spectrum, and which would be the most pessimistic sit- 
uation for producing anisotropy or source correlations in view 
of their large magnetic deflections. 

This paper is laid out as follows. In Section 2 we describe 
our models and methods of calculation and statistical analy- 
sis. Parametrization of photodisintegration cross sections to 
calculate the mean free paths of nuclei in cosmic background 
radiation fields is summarized in Section 2.1. Section 2.2 is 
dedicated to our models for the GMF and EGMF. A method 
of calculation for nuclear propagation is introduced in Section 
2.3, and that for UHECR arrival distribution is presented in 
Section 2.4. The contents of these two subsections are essen- 
tially based on Ref. [62] that dealt with UHE proton propaga- 
tion in intergalactic space, and here we extend their methods to 
heavy nuclei. We describe our statistical methods for studying 
anisotropy and cross correlations with sources in Section 2.5. 
In Section 3, we present our simulation results, compare with 
the current data [34, 63], and discuss the implications thereof 
for the sources of UHE nuclei. We also give special attention to 
Centaurus (Cen) A, for which several PAO events positionally 
correlate [63]. We summarize in Section 4. 

2. Details of Calculations 

2.1. Photodisintegration 

A nucleus emits its constituent nucleons by interactions with 
photons (photodisintegration). Since the number of emitted 
nucleons per interaction is stochastically determined following 
branching ratios dependent on interaction energy, a lot of disin- 
tegration paths are realized during propagation in photon fields. 
However, all the branches in photodisintegration are not well 
known yet. This study adopts a simple disintegration track on 
the nuclear chart proposed by Ref. [26]. There is only one nu- 
clear species for each nuclear mass number A in this treatment. 
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Figure 1 : Energy-loss length of Bethe-Heitler pair creation (red) and mean free 
path of total photodisintegration (green) in a model of cosmic background radi- 
ation field [68] for jgFe- The contribution of individual processes of photodis- 
integration is also shown for comparison. 



Photodisintegration can be phenomenologically classified 
into 4 processes: giant dipole resonance (GDR), quasi-deuteron 
(QD) process, baryonic resonance (BR), and photofragmenta- 
tion (PF). Since the importance of the latter two processes ap- 
pears above 10 22 eV for Fe nuclei in the CMB because of their 
very high threshold energy (fig. 1, and see also Section 2.3), 
our assumption of Ef^ x = 10 215 eV allows to neglect these 
processes. This study adopts the parametrization of photodisin- 
tegration cross-sections developed by Ref. [64]. The mean free 
paths of nuclear photodisintegration in cosmic background pho- 
ton fields can be calculated from the cross-sections (explained 
in Section 2.3). 

GDR has the lowest threshold energy in the photodisintegra- 
tion processes, typically ~ 10 MeV in the nucleus rest frame. 
A nucleus emits several nucleons or a particles, but one nu- 
cleon emission is dominated. The treatment of Ref. [64] ap- 
proximated GDR as a process emitting one nucleon. The cross- 
section of GDR is parametrized by a Lorenzian function 
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for e' > e' A th GDR and otherwise. Here e' is photon energy in 
the nuclear rest frame, &q = 0.72 mbarn, e'gdr = 35.3 MeV, 
and Tgdr = 15.1 MeV. The threshold energy of this process is 
also scaled as a function of A, e' A th GDR = 16.6 A~ 1/6 MeV. 

The quasi-deuteron process is treated as a two nucleon emis- 
sion process. The parametrization of its cross-section is 
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for e' > f' where e' = 33 3 A' 116 MeV 

lOl fc ^ fc A,th,QD W1,C1C fc A,th,QD -->-->-> ^ 1V1CV. 

While the scaling laws on A can well reproduce experimental 
data for large A nuclei, the behavior of small nuclei (A < 10) 



is different because of their small number of nucleons. Thus, 
other analytical functions of the cross-sections are used for 
deuteron, trinucleon (triton or 3 He), helium, and 9 Be instead 
°f o~a,gdr(c'), o~a,qd(€')- ^Be is assumed to decay into two pis 
and one nucleon. We also assume that a nucleus with A = 8 
which is produced by QD of a particle with A = 10 imme- 
diately decay into two !jHe because there is no stable particle 
with A = 8. Therefore, nuclei between A = 5 and A = 8 are not 
observed at the Earth in this treatment. 

The photoabsorption cross-section of a deuteron was derived 
by Ref. [65] as 
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is a photoabsorption cross-section for a zero range nuclear force 
[66]. h, m p , c, a, B d = 2.227 MeV, r eff = 1.79 fm, and 
<TTp = 1-97 x 10~ 4 mbarn are the Planck constant, the proton 
mass, the speed of light, the fine structure constant, the nuclear 
binding energy of a deuteron, the effective range of nuclear 
force, and the cross-section of Thomson scattering for proton 
mass, respectively. After this reaction, a deuteron is broken to 
a proton and neutron. 

A trinucleon has two interaction channels; it emits a nucleon 
and is changed into a deuteron, or three body decay. The cross- 
sections of these processes can be well fitted by using cr BP , 

1.4cr BP (e', 5.8MeV) for T 3 (y, A?)D 2 (6) 
1.7cr BP (e', 7.3MeV) forT 3 (y,3/V) (7) 
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The cross-sections of one or two particle emission of a- 
particle are parametrized as 

<(e0 = (3.6mbarn)Pl(e', 19.8MeV, 27MeV, 5) (8) 

cr^(e') = 1.4cr BP (e', 26.1MeV), (9) 

respectively, where 



Plfc *th, *max, a) = 
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Here, 0(x) is a step function and a — a(x max /x t h). 

The cross-section of photo-absorption of ^Be is parametrized 



as 

where 
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(e) = 1.5<r L (e', 1.7MeV,0.3MeV) 

+ 1.6 <x L (e', lO.OMeV, lO.OMeV) 
+ 3.5o- L (e',25MeV, 15MeV) mbarn (12) 
o- 9 2 Be (e') = 9 cr BP (e', 16.9MeV). (13) 

We should also treat the propagation of protons as secon- 
daries. For protons, photopion production plays an important 
role above 6x 10 19 eV. The cross-section and energy-loss rate of 
protons are calculated by a photomeson production event gen- 
erator SOPHIA [67] . 
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2.2. EGMF and GMF 

Based on the cosmological structure formation theory, den- 
sity fluctuations produced by inflation have gravitationally 
grown and form large-scale structure of the current universe. 
Galaxy surveys have confirmed through galaxy distribution 
that matter distribution is highly structured. Since magnetic 
fields has been also amplified through the compression and 
dynamo mechanism of astrophysical plasmas simultaneously, 
they could be also structured following the matter distribu- 
tion. Several simulations of cosmological structure forma- 
tion with magnetic fields have confirmed this expectation qual- 
itatively [69, 70, 71, 72], but EGMF structures highly de- 
pend on modeling due to few observational constraints. The 
difference appears, for example, in the volume filling factor 
of strong EGMFs (see Ref. [73]). Therefore, a variety of 
EGMF modelings gives us different results on the deflection 
of UHECRs and source identification possibility even for pro- 
tons [69, 70, 71, 74, 75, 76]. In addition to the uncertainty, 
it is technically difficult to calculate the arrival distribution of 
UHECRs with a complex EGMF structure because only a very 
small fraction of UHECRs injected from sources can arrive at 
the Earth. Ref. [74] developed a numerical method to solve 
this problem for UHE protons, but this method cannot apply 
to the cases of heavy nuclei owing to the stochastic nature of 
photodisintegration. In order to avoid the uncertainty of EGMF 
modeling and the numerical difficulty, this study assumes a uni- 
form turbulent EGMF model with the strength of Begmf, the 
correlation length of /Iegmf, and the Kolmogorov power spec- 
trum [62]. The averaged property of EGMFs has been con- 
strained by Faraday rotation measurements of distant quasars 
as S E gmf/Iegmf 1/2 < (10nG)(lMpc) 1/2 [77, 78]. A recent anal- 
ysis of CMB and matter power spectra derives an upper limit 
of the primordial magnetic field at present, which is regarded 
as a magnetic field in voids, of B v < 3 nG (95 % C.L.) for the 
coherent length of A v = 1 Mpc [79]. Several simulations of 
cosmological structure formation with magnetic fields indicate 
B v <sc 1 nG [71, 75]. On the other hand, pair-halo and/or pair- 
echo emission from TeV blazars may provide lower bounds, 
suggesting B v > 10~ 18 - 10~ 17 G from the data of Fermi Large 
Area Telescope [80, 81, 82]. Even magnetic fields in voids 
are negligibly weak, the effective magnetic fields of cosmic 
structures such as filaments and clusters of galaxies may have 
Begmf^egmf 1/2 < 1 nG Mpc 1/2 [83]. In this study B E gmf is 
treated as a parameter, while /Iegmf = 1 Mpc is fixed which is 
motivated by the average separation distance between galaxies. 

The GMF has been relatively well observed compared to 
EGMFs through Faraday rotation measurements of radio pul- 
sars and extragalactic radio sources. It has a regular component 
with a few fiG on average in the Galactic disk, whose shape is 
similar to the spiral structure of the Milky Way. There are also 
turbulent fields with the magnetic strength 0.5-2 times as large 
as that of the spiral component [84]. Its correlation length has 
been lower-limited to Aqmf ~ 100 pc due to the finite resolution 
of radio observations. Both the components significantly con- 
tribute to the total deflections of UHE nuclei [85]. However, at 
present, even the newest data are difficult to discriminate the 
regular structure between the bisymmetric and axisymmetric 



models [86], and the difference produces significantly different 
propagation of UHECRs especially in the case of heavy nuclei, 
e.g., nuclear propagation is sensitive to GMF models [87]. In 
order to avoid the model dependence, we simply assume that 
the deflection angles of UHECRs do not depend on directions. 
In Ref. [88], a GMF model with a weak dipole field predicted 
a typical deflection angle of 40°-50° for irons with ~ 6 x 10 19 
eV on average in the whole sky. Assuming a rigidity scaling 
additionally, we set the typical deflection angles of UHECRs in 
GMF to be 6» G mf = 1.0°Z(£/10 20 eV)- 1 where E is the energy 
of UHECRs. 

2.3. Calculation of Propagation 

The calculation of nuclear propagation can be divided into 
two parts; propagation in extragalactic space and in Galactic 
space. The size of the Galactic space (assumed to be 40 kpc) is 
much smaller than the energy-loss length of Bethe-Heitler pair 
creation. On the other hand, the mean free path of photodis- 
integration is close to the size of Galactic space above 10 205 
eV for jgFe. However, it does not affect the propagation in the 
Galaxy because particles traveled for longer distances than the 
mean free path in extragalactic space before entering the Milky 
Way. Thus, the nuclear reactions can be neglected during the 
propagation in the Galaxy. This allows us to treat the effect 
of the propagation as modifications of the arrival directions of 
UHECRs. The modifications will be described in Section 2.4. 
Below, we introduce a method of the propagation of UHECRs 
in extragalactic space. 

We calculate the propagation trajectories of UHECRs includ- 
ing secondary particles in extragalactic space, and estimate the 
nuclear species, energies, deflection angles of nuclei arriving 
at Earth. In our method, the source of primary nuclei (^Fe) 
is located at the center of coordinates in extragalactic space. 
The primary nuclei are ejected from the source isotropically, 
and then their trajectories are calculated step by step includ- 
ing secondary particles taking magnetic deflections, photodis- 
integration, and Bethe-Heitler pair creation into account. Adia- 
batic energy-loss due to cosmic expansion is neglected because 
the energy-loss rate is much lower than the others. Since this 
study focuses on the highest energy range, E > 5.5 x 10 19 eV, 
the calculation of the trajectories is finished at 500 Mpc from 
the source. We confirmed that nuclei with the energy above 
E > 5.5 x 10 19 eV cannot arrive at the Earth from more distant 
sources even for rectilinear propagation. The nuclear species, 
energy, and deflection angles of the particles are recorded every 
1 Mpc. We divide a logarithmic energy bin into 10 per decade, 
i.e., 10 1975 -10 19 - 85 eV, l0 19 - 85 -10 19 ' 95 eV, • • •, 10 2145 -10 21 - 55 eV, 
and inject 10,000 Fe nuclei in each bin. The record allows us to 
make the distribution of the nuclear species, energy and deflec- 
tion angles of the particles at 500 distances. For instance, when 
we consider UHECRs injected from a source with the distance 
of 100 Mpc, the 100th distributions of them are regarded as 
their distributions at the Earth. 

In the calculation of propagation of nuclei, photo- 
disintegration processes are treated as a stochastic process, i.e., 
by a Monte-Carlo method. The mean free path of photo- 
disintegration process ;' (= GDR, QD) of a nucleus with the 
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nuclear mass number A and the energy E can be calculated by 
using the expressions of cross-sections described in Section 2.1 
as [89], 



Table 1: Local number densities of active objects. 



Objects 



Density [Mpc 3 ] Refs. 



1 



1 f°° de dn y r s — , , 



A A j(E) 8/?£ : 

where (5, e, s, and M A are nuclear velocity normalized by c, 
the energy of photons in the laboratory frame, the Lorentz in- 
variant mass squared, and nucleus mass, respectively. cr A j(s) is 
cross-section of photo-disintegration process considered, s m ; n = 
M A 2 c 4 + 2M A c 2 e' A lhi , s max = M A 2 c 4 + 2Ee(l + 0), and e th = 
(s m ; n - M A 2 c A ) \2E{\ + /?)] _1 . M A is simply defined as M A = 
Am„, where m„ = 931.494 MeV is nucleon mass. dn y jde is the 
number density of extragalactic background light (EBL) pho- 
tons including CMB photons. This study adopts the low-IR 
model of Ref. [68]. Since the Bethe-Heitler pair creation pro- 
cess is treated in a continuous energy-loss approximation (see 
below), the occurrence of reactions is judged by the method 
discussed in Ref. [90]. GDR and QD are treated independently. 
We approximate that the Lorentz factor of nuclei is unchanged 
by photodisintegration. 

For the energy-loss length of the Bethe-Heitler pair creation, 
we adopt an analytical fitting formula developed by Ref. [91], 

j-i= are 2 cZ 2?^ r d /jL(jL\^ (15) 

dt M A J 2 de\2j) k 2 

where r e is the classical electron radius, m e is electron mass, 
and k = e/m e c 2 . <p(k) is an integral similar to the first integral 
of Eq. 14, but is parametrized as 



<P(k) = 
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where c\ 
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0.8048, c 2 = 0.1459, 
-86.07, di = 
78.35, and f 3 



c 3 = 1.137 
50.96, d = 
= 1837. 



k<25 
k>25, 

x 10~ 3 , c 4 
-14.45, d 
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2.4. Calculation of Arrival Distribution 

The arrival distribution of UHECRs can be calculated using 
the distributions obtained by the method written in Section 2.3 
with an assumed source model. 

We adopt a source model developed in Ref. [74], in which 
UHECR sources are distributed following galaxy distribution 
actually observed. The source distribution of this model is 
based on the galaxy catalog of the Infrared Astronomical Satel- 
lite (IRAS) Point Source Redshift Survey (PSCz) [92]. This 
galaxy catalog covers ~ 84% of the whole sky, and is enough 
to construct the model of matter distribution in local Universe 
(e.g., Ref. [71]). Selection effects of galaxies (e.g., unobserved 
sky and galaxies fainter than the flux limit) are compensated by 
a luminosity function [93], and then the model of galaxy distri- 
bution following matter distribution is constructed within 200 
Mpc. Outside 200 Mpc, we assume isotropic distribution due 
to the lack of completeness of galaxies in the catalog. This as- 
sumption does not change the results of this study, since distant 



Seyfert galaxy 
Bright quasar 
Fanaroff-Riley I 
Fanaroff-Riley II 
BL Lac objects 



1 X 10~ 3 
1 x 10~ 6 
8 x 10~ 5 
3 x 10~ 8 
3 x 10" 7 



[100] 
[101] 
[102] 
[103] 
[103] 



sources do not contribute to anisotropy. We select galaxies ran- 
domly from the galaxy distribution following n s , which is the 
number density of UHECR sources as a parameter, and regard 
a set of the selected galaxies as a source distribution. Although, 
n s has been estimated as ~ 10~ 4 Mpc~ 3 [94, 95] on the assump- 
tion of pure proton composition, it is not clear in the case of 
heavy nuclei. A smaller n s may be allowed owing to larger de- 
flections. In this study, we survey over n s = 10 3 -10~ 7 Mpc -3 . 
Corresponding objects are shown in Table 1 . All the analyses 
in this study are limited to steady sources. We assume that all 
the sources emit UHECRs with the same power for simplic- 
ity. The injection spectrum of the particles is assumed to be a 
power-law oc E~ 2X) with E^! dx = 10 21,5 eV. We confirmed that 
this spectrum can reproduce observed spectrum above 6 x 10 19 
eV within uncertainties. 

In Table 1, Fanaroff-Riley (FR) I galaxies and BL Lac objects 
are in different categories. However, following a unification 
theory of radio-loud AGNs, these originate from the same pop- 
ulation with the only difference due to the viewing angle [96]. 
Thus, if the deflection angles of UHECRs during propagation 
in extragalactic space is much larger than a typical opening an- 
gle of their jets (~ 0.1 radian), we can see these objects as the 
essentially same objects in the viewpoint of UHECR sources. 

An assumed n s constrains the energy budget of UHECR 
emission per source. The observed spectrum of UHECRs re- 
quires ~ 10 44 erg Mpc 3 yr 1 [97, 98, 99] for UHECR emission 
above 10 19 eV. Thus, the energetics of UHECR emission re- 
quired per source is ~ 3 x 10 40 (n. s /10 _4 Mpc _3 ) _1 erg s _1 . 

The flux contribution of each source is inversely propor- 
tional to the distance squared. Following the injection spec- 
trum weighted by this factor and the distribution of UHECRs 
in deflections, energies, and nuclear mass number (see Section 
2.3), we can simulate the arrival distribution of UHECRs. In 
addition, we can add the contribution of the GMF in the arrival 
distribution. Once the arrival directions of UHE nuclei are real- 
ized following the method above (considering only EGMF), the 
arrival directions are modified following the two dimensional 
Gaussian distribution with the zero mean and the standard de- 
viation Of #GMF- 

2.5. Statistical Quantities 

Anisotropy in the arrival distribution of UHECRs can be 
quantified by a cumulative auto-correlation function [104], 



CM = 



EE{< 0) - 2EE'(< 0) + E'E'(< 0) 
E'E'(< 0) ' 



(17) 
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Figure 2: Cumulative auto-correlation functions of simulated events for the source models with n s = 10~ 3 , 10~ 4 , 10~ 5 Mpc~ 3 (left) and 10" 6 , 10~ 7 Mpc~ 3 (right). 
The simulated event sets mimic the recently published events of the PAO, i.e., 69 events with the energy above 5.5X10 19 eV, assuming the PAO aperture. Both EGMF 
(Begmf = 1 oG) and the GMF are taken into account. Points are the averaged values of the cumulative auto-correlation functions over 1000 source realizations for 
each n s . The error bars correspond to 68% (lcr) errors. 



where E and E' denote UHECR events and events randomly 
distributed following the aperture of a UHECR detector, respec- 
tively. £■£■(< 9) is the number of UHECR event pairs with the 
separation angle less than 9 divided by the total number of pairs. 
EE'(< 6) and E'E(< 9) are defined similarly to EE{< 9). This 
function is essentially equivalent to the quantity that the num- 
ber of UHECR event pairs with the separation angle within 9 
is divided by the total number of pairs and by a solid angle 
within 9 if the aperture of UHECR observatories is uniform. 
However, since, in general, the aperture is not uniform in the 
cases of ground-based detectors, the non-uniformity should be 
corrected. E' can correct this in this formula. This function 
was originally proposed by Ref. [105] as a differential auto- 
correlation function. This correction simplifies the interpre- 
tation of the cumulative auto-correlation function; C e {9) > 
shows the positive excess of UHECR distribution within 9 com- 
pared to isotropic distribution. 

Many tests of the correlation between the arrival directions of 
UHECRs and the positions of their source candidates are based 
on the excess of the number of events within a circle centered 
by the source candidates over randomly distributed events. This 
study also adopts a similar method to study the correlation, us- 
ing a cumulative cross-correlation function, which is defined 
similarly to the cumulative auto-correlation function, 

ES (< 9) - ES '« 9) - E'S (< 8) + E'S '(< 9) 
' E'S'(<8) 

where S and S ' denote UHECR sources and sources randomly 
distributed following selection effects of the catalog of the 
source candidates. S' is isotropically distributed in the entire 
sky in this study because of the correction of selection effects 
(see Section 2.4). ES(< 9) is the normalized number of pairs 
between UHECR events and sources with the separation an- 
gle less than 9. ES'(< 9), E'S(< 9), and E'S'(< 9) are simi- 
larly defined. The interpretation of this quantity is also simple; 
C es {9) > means positive correlation within 9. 



In order to calculate Eqs.17 and 18, the distribution of E' 
should be understood. The aperture of a ground array depends 
on the declination of observed directions reflecting the daily ro- 
tation of the Earth. The declination dependence of the aperture 
can be analytically estimated as [106] 

lo(6) cc cos(flo) cos(<5) sin(a m ) + a m sin(an) sin((5), (19) 



where a m is given by 



and 



if ^ > 1 
if £ < -1 

otherwise 



cos(©) - sin(an) sin((5) 




cos(ao) cos(<5) 



(20) 



(21) 



when observation time is sufficiently larger than a day. Here, ao 
is the celestial latitude of the ground array and is the zenith 
angle for a data quality cut because of experimental reasons. 
For the PAO, a = -35.2° and = 60° [33, 34]. The E' sample 
is generated following this equation. We set the number of E' 
events to be 200000 in order that the distribution of random 
events reflects the PAO aperture sufficiently. Also, we checked 
that this particular choice of the number of random events does 
not affect results. 

3. Results 

3.1. Anisotropy 

Fig. 2 shows the cumulative auto-correlation functions of 
simulated events with 68% errors in the cases of n s = 10 -3 - 
10" 7 Mpc~ 3 . They are divided into two panels for visibility. 
The simulated event set imitates the recent PAO data [107], 
e.g., 69 events above 5.5 x 10 19 eV are distributed following 
the aperture geometry of the PAO. Begmf = 1 nG and GMF are 
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Figure 3: Cumulative auto-correlation functions of 69 events above 5.5 X 10 
eV on the assumption of the PAO aperture for n s = 10~ 4 Mpc~ 3 and Begmf = 1 
nG. The GMF is considered (green) and is not considered (red) for comparison. 



Figure 4: Cumulative auto-correlation functions of 69 events above 5.5 X 10 
eV simulated on the assumption of the PAO aperture for different Begmf, 0.0, 
0.1, and 1 nG for comparison. n s = 1CT 5 Mpc~ 3 and the GMF are considered. 



taken into account. The points and error bars of these cumula- 
tive auto-correlation functions are calculated as follows. First 
of all, 69 events are generated from a source distribution with 
the source number density of n s , and then a cumulative auto- 
correlation function is calculated from this event set. Next, the 
former step is repeated over 1000 source distributions with the 
same n s . Then, the averages and errors are calculated from the 
1000 cumulative auto-correlation functions. When the number 
of events is small, the distribution of the values of the auto- 
correlation functions in a bin are asymmetric and therefore the 
standard deviation of the 1000 auto-correlation functions is not 
a good indicator of errors. In order to estimate anisotropic er- 
ror bars, we arrange the values of the 1000 cumulative auto- 
correlation functions in each bin in order, calculating the range 
in which the central 68% of the values are included, and then 
this range is regarded as a 68% (lcr) error. This error includes 
not only an error due to the finite number of events (Poisson er- 
ror) but also an error due to the uncertainty originating from the 
selection of UHECR sources. Finally, the averages and 68% er- 
rors are plotted in the figure. All the figures in this paper follow 
this rule, as long as there is not exceptional description. 

All the auto-correlation functions except for n s = 10~ 7 
Mpc 3 are consistent with zero within 68% errors in this an- 
gular range, i.e., the arrival distributions are consistent with 
isotropic distribution. Note that the case of n s = 10~ 6 Mpc 3 
corresponds to only 4 sources within nearby 100 Mpc. In spite 
of such a small number, the predicted arrival distribution is 
consistent with isotropy because of the strong effects of the 
GMF and EGMF These magnetic fields make even the case of 
n s = 10~ 7 Mpc~ 3 become consistent with isotropy within 95% 
error at small angular scale (< 20°). 

The GMF significantly contributes to the arrival distribution 
of UHECRs. Fig. 3 shows the cumulative auto-correlation 
functions of 69 simulated events in the case of n s = 10 Mpc 3 
and Begmf = 1 nG. In order to see the effect of the GMF, one 
(green) takes the GMF into account, but the other (red) does 
not consider the GMF. Only in this figure, the medians of 1000 



cumulative auto-correlation functions are plotted instead of the 
averages for visibility because the very large values of cumula- 
tive auto-correlation functions are realized in some realizations 
by very nearby sources at small angular scale and the averaged 
values of the auto-correlation functions become larger than the 
upper edge of the 68% error bars when the GMF is not taken 
into account. Note that the medians will be close to the av- 
erages when the number of events increases. Focusing on the 
first bin, we can see that the medians in both the cases are close 
to the lower edges of the error bars. This means that source 
distributions including nearby sources accidentally make strong 
anisotropy whereas no event pair appears in the first bin in many 
realizations because of the small solid angle within 1° and the 
small number of simulated events. Thus, the distribution of the 
values of the cumulative correlation functions in the first bin has 
a long tail to large values. Although both the cumulative auto- 
correlation functions are consistent each other (and also with 
isotropic distribution) within 68% errors, the error bars are sup- 
pressed when the GMF is considered, that is, the probability 
that large values of the auto-correlation function are realized 
becomes smaller. The GMF strongly weakens anisotropy pro- 
duced by nearby sources. 

Furthermore, we can find that the GMF can dominantly con- 
tribute to suppressing anisotropy. Fig. 4 shows the dependence 
of the cumulative auto-correlation functions on Begmf in the 
case of n s = 10~ 5 Mpc~ 3 and considering the GMF. For com- 
parison, Begmf = 0, 0.1, and 1 nG are adopted. The three cu- 
mulative auto-correlation functions are quite similar, and there- 
fore are controlled by the GMF. This situation is the same down 
to the case of n s = 10~ 6 Mpc -3 . Small difference of the aver- 
ages among the cases of the three EGMF strengths appears only 
for /I., = 10~ 7 Mpc -3 ; Begmf 0.1 nG predicts a bit higher val- 
ues of the averages, but the difference is still within the error 
bars. Thus, the effect of the GMF is essential and should be 
inevitably considered when the anisotropy of UHECR arrival 
distribution is tested. 

Here, we are interested in how often a positive value of the 
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Figure 5: Probability that the positive values of the cumulative auto-correlation functions are not realized for 69 events above 5.5 X 10 eV simulated on the 
assumption of the PAO aperture. The upper two panels and the lower two panels assume Begmf = 0.0 nG and 1.0 nG, respectively. Only the right two panels take 
the GMF into account. The probability that the positive excess of events appears is indicated in the panels. 
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auto-correlation function is realized in each bin. This provides 
the probability that the positive excess of events appears in each 
angular scale. This probability is calculated as follows. We 
simulate a set of UHECR events from a source distribution and 
calculate the cumulative auto-correlation function of them. This 
calculation is repeated over 1000 source distributions with the 
same n s . Then, we focus on the values of the 1000 cumulative 
auto-correlation functions in an angular bin, counting the num- 
ber of the auto-correlation functions which are positive, and the 
number is divided by 1000. This step is repeated over all the an- 
gular bins. Since only 1000 source distributions are considered 
(because of the computational limitation), it cannot resolve the 
probability less than 10~ 3 , but it is enough for this study. By 
definition, isotropic distribution predicts this value of ~ 50% 
due to statistical fluctuation. This probability also can be re- 
garded as a confidence level in the one-side statistical test to 
rule out isotropic distribution. For visibility, we plot one minus 
this probability throughout the paper, which is the probability 
that the positive excess of events is not realized. 

Fig. 5 shows such probabilities for 69 simulated events above 
5.5 x 10 19 eV. The aperture geometry of the PAO is applied. The 
upper two panels and the lower two panels are calculated under 
the conditions of Begmf = nG and 1 .0 nG, respectively. The 
GMF is considered for only the right two panels. 

In the case of no magnetic field (upper left), the strong posi- 
tive excess of events appears at small angular scale. This is nat- 
ural because the propagation trajectories are not deflected. In 
general, the probability that the positive excess of events is re- 
alized in a certain bin is higher for a smaller n s at small to inter- 
mediate angular scale because each source has stronger power 
and stronger clusterings of UHECR events are expected in the 
directions of their sources. The positive anisotropy is produced 
at more than 99% for n s < 10~ 5 Mpc~ 3 at small angular scale. 
Note that this does not mean that the positive excess signal of 
a set of observed events compared to isotropic distribution is 
larger than 99% confidence level at the angular scale. 

Magnetic fields change the shape of the probability curves. 
When the GMF is taken into account (upper right), the angu- 
lar scale at which the value of a probability curve is minimized 
moves from the first angular bin to larger angular scale for all 
the n s . The possibility of positive excess is maximized at this 
angular scale. This angular scale reflects the deflection angles 
of UHECR trajectories by the GMF. Since these are indepen- 
dent of the distances of sources, it does not depend on n s . Al- 
though the EGMF gives a similar effect to the GMF (lower left), 
the angular scale of the minimum of a probability curve is dif- 
ferent and depends on n s because the deflections of UHECRs by 
the EGMF are dependent on their propagation distance and the 
nearest sources, which contribute to anisotropy most strongly, 
are more distant in the case of a smaller n s . In both cases only 
sources with n s = 10~ 7 Mpc 3 produce positive excess in inter- 
mediate angular scale with more than 95 % probability. 

The result when both the GMF and EGMF are taken into 
account (lower right) can be interpreted as combination of the 
upper right and lower left panels. The contribution of the two 
magnetic fields reduces the probability that positive excess of 
events appears and moves the angular scale of the maximum of 
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Figure 6: Same as the lower right panel of Fig. 5, but for 1000 events above 
5.5 X 10 19 eV simulated on the assumption of a uniform aperture. 



the probability to larger angular scale. 

At the end of this section we consider the case with the num- 
ber of events expected by future UHECR experiments such as 
the Northern Pierre Auger Observatory [108] and Extreme Uni- 
verse Space Observatory on board Japanese Experiment Mod- 
ule (JEM-EUSO) [109]. Fig. 6 is the same figure as that in the 
lower right panel of Fig. 5 (the GMF and EGMF are consid- 
ered), but for 1000 events above 5.5 x 10 19 eV on the assump- 
tion of the uniform aperture. All the cases of n s predict the 
positive excess of events at intermediate angular scale at more 
than 99 %. We notice that the case of n s = 10~ 3 Mpc~ 3 shows 
a larger probability of the positive excess than that of n s = 10~ 4 
Mpc~ 4 at around 40°. The extremely close sources included 
for n s = 10~ 3 Mpc~ 3 (within a few Mpc) causes this feature. 
When sources within 5 Mpc are artificially neglected for exam- 
ple to check this, the probability curves for n s > 10" 4 Mpc~ 3 
change significantly and the inversion disappears. In this case, 
the probability for n s = 10~ 3 Mpc~ 3 becomes a bit smaller to 
90%. 

3.2. Correlation with sources 

Fig. 7 shows the cumulative cross-correlation functions be- 
tween the arrival directions of 69 UHECRs above 5.5 x 10 19 eV 
simulated on the assumption of the PAO aperture and the celes- 
tial positions of their sources within 75 Mpc with 68% errors in 
the cases of n s = 10 _3 -10 5 Mpc -3 . The 75 Mpc is motivated by 
the first report of the correlation between UHECRs and nearby 
extragalactic sources by the PAO [33]. Begmf =1-0 nG and 
the GMF are taken into account. The error bars and the aver- 
aged points of the cross-correlation functions are estimated sim- 
ilarly to those of cumulative auto-correlation functions. Note 
that we do not plot the cumulative cross-correlation functions 
for n s = 10~ 6 and 10~ 7 Mpc~ 3 because there are source dis- 
tributions not having sources within 75 Mpc in 1000 source 
realizations due to the small number densities. The cumula- 
tive cross-correlation functions for n s = 10~ 3 and 10~ 4 Mpc~ 3 
are fully consistent with zero within 68% errors in this angu- 
lar scale, i.e., are consistent with no correlation. The case of 
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Figure 7: Cumulative cross-correlation functions between the arrival directions 
of simulated UHECR events and the celestial positions of their sources with 
n s = 10~ 3 , 10 -4 , and 10~ 5 Mpc~ 3 within 75 Mpc. The simulated event sets 
imitate the published events of the PAO, i.e., 69 events with energies above 
5.5 X 10" eV simulated on the assumption of the PAO aperture. Both EGMF 
(Segmf = 1 nG) and the GMF are taken into account. The points are the 
averaged values of the cumulative cross-correlation functions over 1000 source 
realizations for each n s . The error bars correspond to 68% errors. 



n s - 10~ 5 Mpc~ 3 is positive at small angular scale with 68% 
errors, but is consistent with no correlation within 95% errors. 
The error bars are longer for smaller n s because the fluctuation 
of source positions is larger. 

The cross-correlation between UHECRs and their sources is 
also dominantly affected by the GMF. Fig. 8 is similar to Fig. 
4, but the cumulative cross-correlation functions between the 
69 UHECRs and their sources with n s = 10~ 4 Mpc 3 within 
75 Mpc. Begmf = (red), 0.1 (green), and 1 nG (blue) are 
considered with the GMF for comparison. Similarly to Fig. 4, 
the three cumulative cross-correlation functions are very sim- 
ilar. The situation is the same even for n s = 10 3 and 10~ 5 
Mpc 3 . Even in this figure, we can notice the importance of the 
GMF on the discussions of the arrival distribution of UHECRs. 

Estimating the probability that positive correlation between 
UHECRs and their sources appears is also useful. We can cal- 
culate this probability by a similar method to calculate Fig. 5. 
Following that figure, we plot the possibility that the values 
of the cross-correlation function are not positive for visibility 
in Fig. 9 similarly to Fig. 5. We consider 69 UHECRs above 
5.5 x 10 19 eV on the assumption of the PAO aperture. The 
sources used for calculating the cross-correlation function are 
within 75 Mpc from the Galaxy. The basic tendency of the 
probability curves is similar to Fig. 5. In the case of no mag- 
netic field, strong correlation is predicted at the smallest angu- 
lar bin because of the absence of the deflections of UHECRs. 
Taking the GMF and/or EGMF into account, the probability 
curves have a minimum at intermediate angular scale. For the 
cross-correlation function, the angular scale at which gives the 
minimum depends not only the deflection angles of UHECRs 
but also the typical separation angle between sources. The po- 
tential minimum is at ~ 20° for n s = 10~ 5 Mpc~ 3 if the GMF 
and/or EGMF are considered, but the probability of the positive 



Figure 8: Cumulative cross-correlation functions between 69 UHECRs above 
5.5 X 10" eV simulated under the PAO aperture and sources within 75 Mpc for 
different Begmf, 0.0, 0.1, and 1.0 nG for comparison. «, = 10~ 4 Mpc~ 3 and 
the GMF are considered. 



correlation is less than 95%. 

Again, we demonstrate the prospects of future UHECR ob- 
servatories. Fig. 10 shows cumulative cross-correlation func- 
tions (left) and the corresponding probability curves (right). We 
consider 1000 UHECRs above 5.5 x 10 19 eV on the assumption 
of the uniform aperture and both the GMF and EGMF (Begmf = 
1 nG). In the left panel, the cumulative cross-correlation func- 
tions in the cases of n s - 10~ 4 and 10~ 5 Mpc~ 3 are clearly sepa- 
rated from zero within 68% errors, in comparison to Fig. 7. Al- 
though it is difficult to see, even the cases of n s = 10~ 3 Mpc~ 3 
is also inconsistent with no correlation within 68% errors. The 
probability that the positive correlation appears is shown in the 
right panel. The cases of n s — 10~ 4 and 10~ 5 Mpc~ 3 produce the 
correlation with more than 99 % at small angular scale. Since 
stronger correlation is expected for the smaller n„ next gen- 
eration UHECR experiments can observe positive correlation 
between UHECRs and their sources at small angular scale with 
more than 99 % probability if n s < 10~ 4 Mpc~ 3 . 

3.3. Cen A as Nearest UHECR Source 

The recent PAO data with energies above 5.5 x 10 19 eV 
show the overdensity of the arrival directions in the direc- 
tion of Cen A and that 2 events positionally correlate with 
the nucleus position of Cen A [107]. On the other hand, Cen 
A has been expected as nearby production site of UHECRs 
[110, 111, 112, 113, 114]. In this section, we specially focus 
on Cen A, and discuss the possibility that Cen A is the nearest 
UHECR source from the viewpoint of the arrival distribution of 
UHECRs. 

For this purpose, we use two source models which are 
slightly modified versions of the source model introduced in 
Section 2.4. In the first model, we artificially neglect sources 
within 5 Mpc from the original source model and add Cen A. 
In this source model, all the source distributions include Cen 
A as the nearest UHECR source. We assume that the positions 
and distance of Cen A are (I, b) = (319°, 19.5°) and 4.1 Mpc 
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Figure 9: Probability that the values of the cumulative cross-correlation function between the arrival directions of 69 UHECRs above 5.5 X 10 19 eV simulated under 
the PAO aperture and the celestial positions of their sources within 75 Mpc is not positive. The upper two panels and the lower two panels assume Begmf = 0.0 nG 
and 1.0 nG, respectively. Only the right two panels take the GMF into account. The probability that positive correlation appears is indicated in the panels. 
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Figure 10: (left): Cumulative cross-correlation functions between 1000 UHECRs above 5.5 X 10 eV simulated on the assumption of a uniform aperture and the 
celestial positions of their sources within 75 Mpc. EGMF (Begmf = 1.0 nG) and the GMF are taken into account, (right): Probability that the positive values of the 
cumulative cross-correlation functions shown are not realized. The probability that positive correlation appears is indicated in the panel. 



11 



69 events 




Figure 11: Arrival distribution of UHECRs realized from a specific source distribution including Cen A with n s = 10~ 4 Mpc~ 3 . 69 events above 5.5 X 10" eV 
under the PAO aperture (upper left) and 1000 events above 5.5 X 1() 19 eV under a uniform aperture (lower right) are considered. £^ e ax = 10 21 5 eV is assumed. The 
composition of the UHECRs is represented in color. The source distribution used in this calculation is shown in the lower right panel, but only the sources within 
75 Mpc are plotted as circles for visibility. The radii of the circles are inversely proportional to the distances of sources. The position of Cen A is the center of the 
largest circle. 



following parameters listed in the IRAS PSCz catalog [92]. The 
second source model is for reference, i.e., we artificially neglect 
sources within 5 Mpc from the original model but do not add 
Cen A. Comparing the results from these two source models, 
we can investigate the effect of Cen A to UHECR arrival distri- 
bution. 

Fig. 1 1 demonstrates the arrival distributions of UHECRs 
above 5.5 x 10 19 eV from a specific source distribution includ- 
ing Cen A with n s = 10~ 4 Mpc~ 3 (lower right). This number 
density is comparable with the local number density of FR I 
galaxies, which Cen A is classified into [102]. The source dis- 
tribution is shown as circles in the Galactic coordinates. The 
radii of the circles are inversely proportional to the distances of 
sources. Although we use all the sources up to 500 Mpc in the 
calculations, only sources within 75 Mpc are plotted as circles 
for visibility. The position of Cen A can be identified as the cen- 
ter of the largest circle. The simulated events are 69 under the 
PAO aperture {upper left) and 1000 under the uniform aperture 
(lower left). The colors represent the composition of UHECRs; 
protons (black), helium (blue), light nuclei with 3 < Z < 7 
(light blue), intermediate mass nuclei with 8 < Z < 20 (green), 
and heavy nuclei with 21 < Z < 26 (red), where Z is the nuclear 
number. 

Even the 69 events remind you that a significant fraction of 
the whole events is distributed around Cen A. In this case, about 
30% of UHECRs originates from Cen A. Nevertheless, the ar- 
rival directions of the UHECRs are scattered by GMF (see also 
Ref. [85]). Several black points are secondary protons. Since 
protons are minimally affected by the GMF, they can point out 



their sources. Three protons correlate with the position of Cen 
A very well. If Cen A is a true source, it is possible that the two 
PAO events correlating with the position of Cen A are protons, 
because several plausible GMF models predict very small de- 
flection angles (~ 1°) in the direction of Cen A [115]. Note that 
secondary protons above 5.5 x 10 19 eV require the maximum 
acceleration energy of irons at sources above 10 215 eV. 

The reproducibility of the auto-correlation function of the ob- 
served data by the source model including Cen A should be 
checked. Fig. 12 shows the cumulative auto-correlation func- 
tions of the arrival directions of 69 UHECRs above 5.5xl0 19 eV 
simulated on the assumption of the PAO aperture. The source 
models with n s = 10~ 4 Mpc 3 explained above are applied. The 
error bars represent 68% (thick) and 95% (thin) errors. Both 
the GMF and EGMF (Begmf = 1 nG) are taken into account. 
The cumulative auto-correlation function of the PAO events is 
also shown for reference. This figure shows that the predictions 
of both source models are consistent with the recent PAO data 
within 95% errors, even in the case that Cen A is the nearest 
UHECR sources. On the other hand, the arrival distribution of 
the recent PAO data is also consistent with isotropic distribu- 
tion above 5.5 x 10 19 eV at 95% confidence level [63]. Thus, 
both possibilities, Cen A is a source or isotropic, are allowed 
at present. Note that the meaning of the 95% in the PAO anal- 
ysis is different from that of the 95% in this analysis. The ob- 
served arrival distribution is consistent with the predictions of 
our source model at 95% confidence level in this study, while 
the distribution is consistent with isotropic distribution at 95% 
confidence level in the analysis of the PAO. 
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Figure 12: Cumulative auto-correlation functions of the arrival directions of 
69 UHECRs with the energy above 5.5 X 10 19 eV simulated under the PAO 
aperture. The source models including (red) and neglecting (green) Cen A with 
n s = 10~ 4 Mpc~ 3 are used. Both EGMF (B EG mf = 1 nG) and GMF are taken 
into account. The error bars represent 68% (thick) and 95% (thin) errors. The 
cumulative auto-correlation function calculated from the recent published PAO 
data [107] is also shown for comparison, (black). 



The PAO data also reported significantly positive excess of 
UHECRs in the arrival distribution around the nucleus position 
of Cen A [107]. This excess gives information to test whether 
Cen A significantly contributes to the total UHECR flux or not. 
Fig. 13 shows the normalized cumulative numbers of events 
used in Fig. 12. The average of the cumulative numbers is 
represented as blue lines with 68% (green) and 95% {light blue) 
errors. The event distribution of the 69 PAO events (black) and 
that of isotropic events (red) are also shown for comparison. 

The source model without Cen A can well reproduce the 
event distribution of the PAO events above 30° within 95% error 
bars (left). The prediction is also consistent with the distribution 
of isotropic events. However, this source model cannot repro- 
duce the excess around Cen A despite that our source model 
includes the overdensity of the Centaurus region. On the other 
hand, the source model including Cen A reproduces the excess 
around Cen A (right). Although the predicted distribution can- 
not reproduce the observed event distribution at large angular 
scales, it indicates that the contribution of Cen A is important 
to reproduce this excess. 

In order to reproduce this observed distribution over all the 
angular scale, the contribution of Cen A has to be reduced or the 
deficit of events should exist at large angular scale from Cen A. 
Possible solutions of the former idea are (i) UHECR flux of Cen 
A is weaker than that of the other FR I galaxies, (ii) there is an 
additional source population with n s > 1CT 4 Mpc~ 3 like Seyfert 
galaxies [7] and the relative contribution of Cen A is reduced, 
(iii) GMF can work more efficiently than what we consider in 
this study. For the latter idea a possibility is Cen A mainly 
contributes to UHECRs in the whole sky (at least southern sky) 
but the UHECRs are not completely isotropized with keeping 
consistency with the observed isotropy mentioned above. In 
any case, scans over large parameter space are needed, which 
is beyond the scope of this paper. Here, we only claim that the 



absence of Cen A cannot reproduce the observed excess. 

The increasing number of events reveals anisotropy by 
nearby sources more clearly, as seen in the lower left panel of 
Fig. 11. If Cen A significantly contributes to the total flux 
of UHECRs, anisotropy in the arrival distribution at angular 
scale of < 30° clearly appears and the order of tens secondary 
protons correlates with the nucleus position of Cen A. In this 
case, ~ 16% of the events originates from Cen A. In addition, 
a significant fraction of secondary protons positionally corre- 
lates with nearby sources other than Cen A. These secondary 
protons are an indicator of nearby sources even if only irons are 
accelerated at sources. 

Fig. 14 shows the cumulative auto-correlation functions 
(left) and cumulative cross-correlation functions (right) of 1000 
events above 5.5 x 10 19 eV simulated on the assumption of the 
uniform aperture. The source models including Cen A (red) 
and neglecting Cen A (green) with n s = 10~ 4 Mpc ~ 3 are ap- 
plied, and both EGMF (Begmf = 1 nG) and GMF are taken into 
account. The two source models are well separated in the cu- 
mulative auto-correlation functions, i.e., anisotropy in events 
detected in the future will be a good clue to understand the 
contribution of nearby UHECR sources. On the other hand, 
although the two source models are not separated in the cu- 
mulative cross-correlation functions as good as in the auto- 
correlation functions, both the cases predict significantly posi- 
tive correlation between UHECRs and their sources at small an- 
gular scale, clearly. Even if primary UHECRs are purely irons, 
information on the location of nearby sources can be extracted 
from the future UHECR data. 

4. Discussion & Summary 

We have investigated the anisotropy in the arrival directions 
of UHECR nuclei and their positional correlations with source 
objects, on the basis of propagation calculations conducted un- 
der various assumptions for the source properties and cosmic 
magnetic fields. Although E^ x = 10 21 5 eV was fiducially 
assumed in the previous section, the results depend somewhat 
on E^ x . Since only iron nuclei above 10 2L5 eV can produce 
secondary protons above 5.5 x 10 19 eV ( fi^ = = 
5.6 x 10 19 [ E^JIO 215 eV ] eV ), the threshold energy for the 
published PAO events and also adopted in this study, secondary 
protons do not contribute to the total flux of UHECRs above it 
if E^ x < 10 21 - 5 eV. Fig. 15 shows the probability curves as 
discussed in Fig. 6, but for E^ x = 10 20 5 , 10 210 , and 10 215 eV, 
in order to check the dependence on £^ e ax . The probability for 
the appearance of a positive excess of events is the highest for 
Zirnax = 10 215 eV at small angular scales, and the angular scale 
of the peak is smaller for higher Ef^ dx . These reflect the fact that 
secondary protons contribute only in the case of E^ ax = 10 21 ' 5 
eV. Furthermore, the probability is slightly higher at small an- 
gular scales for E^ ax = 10 210 eV than for E^ x = 10 20 5 eV 
because light and/or high energy nuclei are included above the 
threshold in the former case. Therefore, positive anisotropy at 
small angular scales is weakened when E^ ax < 10 21 5 eV, owing 
to the lack of secondary protons above threshold, but the prob- 
ability at intermediate angular scales is still greater than 95%. 
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Figure 13: Normalized cumulative numbers of 69 events with energies above 5.5 X 10 eV simulated under the PAO aperture as a function of angular distance from 
the nucleus position of Cen A. In simulations, source models with n s = 1CT 4 Mpc~ 3 including Cen A (right) and neglecting Cen A (left) are used. Both EGMF 
(Begmf = 1 nG) and the GMF are taken into account. The blue lines are the averages of the cumulative numbers of the simulated events. The green and light blue 
regions represent 68% and 95% errors, respectively. For comparison, the event distribution of the 69 PAO events listed in [107] (black) and that of isotropic events 
(red) are also shown. 



We have considered a rather extreme case where only 
iron nuclei are emitted from sources. However, other nuclei 
and/or protons may also plausibly constitute the composition 
of UHECR sources (e.g., a mixed composition model [116] or 
a Wolf-Rayet star model [117]). The other elements can af- 
fect the spectral shape as mentioned in Section 1 due to dif- 
ferent energy losses and magnetic deflections. Although the 
source composition is highly uncertain at the moment, further 
studies of anisotropy and correlations under more realistic as- 
sumptions for the source composition are in order. For instance, 
assuming a mixed composition and a simple rigidity scaling for 
the maximum energies, protons can produce anisotropy at low 
energies comparable to or stronger than the anisotropy due to 
heavy nuclei at the highest energies, because a proton with en- 
ergies Z times smaller than an iron nucleus propagates along 
the same trajectory [118]. Thus, if Cen A is an UHECR source, 
anisotropy at low energies may also be expected. However, 
PAO has reported the lack of significant anisotropy at low en- 
ergies toward the direction of the excess highest energy events 
[119]. This result implies that 1) a significant fraction of UHE- 
CRs may be nuclei even at lower energies, 2) the anisotropy 
at the highest energies is actually weaker than reported, 3) the 
composition at high energies is not as heavy nuclei-dominant as 
indicated from the Xmax measurements, or some combination 
of these three possibilities. Furthermore, strong magnetic fields 
could modify the theoretical expectation through the extension 
of the propagation distance as mentioned in [119], so magnetic 
fields may also be important in this context. 

Although we have adopted simple models for the GMF and 
EGMF to avoid the uncertainties with their structure, the details 
of these fields can actually play important roles in the propaga- 
tion of UHECRs and their arrival directions at Earth, as has 
been well studied for UHE protons [69, 70, 71, 74, 75, 76]. In 
particular, as discussed in Ref. [71, 74, 76] for protons, the 
structures actually observed near the Galaxy are also relevant 



when searching for the evidence of nearby UHECR sources. 
Whether the positional information in the sky is lost or not dur- 
ing propagation depends on the EGMF strength and structures 
in local Universe. The effect of the GMF is inevitable, but the 
deflection angles and directions of UHECRs strongly depend 
on the region of the sky, reflecting the structure of the GMF, 
as studied in Refs. [88, 85] for iron nuclei. Thus more realistic 
GMF models should be applied for correlation studies, although 
such GMF models also have significant uncertainties. The de- 
pendence of the anisotropy signals on different models for the 
GMF and EGMF will be investigated in future work. 

In addition to Cen A, M87 is also an interesting object as an 
UHECR source candidate. Since it is more distant than Cen 
A (~ 16 Mpc), the contribution of M87 to the anisotropy is 
smaller, depending on the EGMF strength. If the power of 
UHECR injection from M87 is comparable to that of typical 
FR I galaxies, it is not expected to contribute strongly to the 
anisotropy. However, if its power is an order of magnitude 
stronger or n s is smaller, anisotropy can appear in the direc- 
tion of M87 even for Begmf = 1 nG. Secondary protons should 
also be a good source indicator in this case. 

This study has assumed the sources of UHECRs to be steady. 
On the other hand, it is also very much possible that transient 
phenomena like GRBs or AGN flares produce UHECRs. Tran- 
sient scenarios are especially relevant for accelerating protons 
up to 10 20 eV in outflows because not many steady sources 
satisfy the luminosity requirement for UHECR acceleration 
[5, 120]. For transient sources, the arrival distribution of UHE- 
CRs will be different from that in the case of steady source sce- 
narios. The time delay of UHECRs by the GMF and EGMFs 
leads to apparent durations of UHECR bursts longer than their 
intrinsic duration. If the apparent duration is longer than the 
observation timescale of human beings, we would misperceive 
an UHECR burst as a steady source [121, 99, 83]. However, 
since the arrival time of UHECRs depends on the UHECR ener- 
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Figure 14: (left): Cumulative auto-correlation functions of the arrival directions of 1000 simulated events above 5.5 X 10" eV. (right): Cumulative cross-correlation 
functions between 1000 simulated events above 5.5 X 10 19 eV and their sources within 75 Mpc. In both panels, a uniform aperture, EGMF (Begmf = 1 nG) and the 
GMF are assumed. The error bars represent 68% errors. 
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gies, features of anisotropy may be strongly energy-dependent 
[121, 83]. This will be evidence for the transient generation 
of UHECRs. Even in transient scenarios, anisotropy of (sec- 
ondary) protons at small angular scale should still be a strong 
hint of UHECR sources. 

To summarize, we have calculated the propagation of UHE 
nuclei above 5.5 x 10 19 eV, taking into account magnetic fields 
in the universe on the assumption that only iron nuclei are in- 
jected by steady sources, and examined the anisotropy in the 
distribution of UHECR arrival directions and positional corre- 
lation between UHECRs and their sources in local universe. We 
found that both the anisotropy and correlation are dominantly 
suppressed by the GMF rather than by EGMFs. Assuming the 
current PAO status (69 events above 5.5 x 10 19 eV), the arrival 
distribution of UHECRs is consistent with isotropy within 95% 
errors for n s > 10~ 6 Mpc -3 and is consistent with no correla- 
tion within 95% errors for n s > 1(T 5 Mpc 3 when the GMF 



and EGMF are taken into account. We also discussed whether 
future experiments improve these situation. 1000 events above 
5.5 x 10 19 eV in the whole sky would reveal positive anisotropy 
with more than 99% probability at intermediate angular scales 
even for n s ~ 10~ 3 Mpc -3 , and correlation between UHECRs 
and their sources in local universe would appear with more than 
99% at small angular scales for n s < 10~ 4 Mpc -3 . In addition, 
we found that the contribution of Cen A is required to repro- 
duce the observed UHECR excess around the Centaurus region. 
Secondary protons from primary heavy nuclei are working pos- 
itively to produce anisotropy at small angular scales, which 
would provide a strong hint of the source location, depending 
on the maximum energy at the UHECR source. 
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